home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_8 / a8_4.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  3.6 KB  |  146 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 8.4, (Steepest Descent  or Gradient Method).
  9. % Section    8.1,    Minimization of a Function, Page 418
  10. echo on; clc; format short; hold off; clear
  11.  
  12. % This program finds the minimum of a function Fn(X),
  13.  
  14. % where X is an n-dimensional vector.
  15.  
  16. % The technique is the gradient method for n-dimensions.
  17.  
  18. % The function  Fn(X)  is to be placed in two M-files.
  19.  
  20. % The function f.m is for graphing and Fn.m is for computing.
  21.  
  22. % The gradient of Fn(X) is to be placed in the M-file Gn.m
  23.  
  24. pause % Press any key to continue.
  25.  
  26. clc;
  27. % function z = f(x,y)
  28. % z = x.^2 - 4.*x + y.^2 - y - x.*y;
  29.  
  30. % function z = Fn(X)
  31. % x = X(1); y = X(2);
  32. % z = x^2 - 4*x + y^2 - y - x*y;
  33.  
  34. % function Z = Gn(X)
  35. % x = X(1); y = X(2);
  36. % G = [2*x-4-y, 2*y-1-x];
  37. % if norm(G)~=0,
  38. %   Z = -G/norm(G);
  39. % else
  40. %   Z = -G;
  41. % end;
  42.  
  43. delete f.m
  44. diary f.m; disp('function z = f(x,y)');...
  45.            disp('z = x.^2 - 4.*x + y.^2 - y - x.*y;');...
  46. diary off;
  47.  
  48. delete Fn.m
  49. diary Fn.m; disp('function z = Fn(X)');...
  50.             disp('x = X(1); y = X(2);');...
  51.             disp('z = x^2 - 4*x + y^2 - y - x*y;');...
  52. diary off;
  53.  
  54. delete Gn.m
  55. diary Gn.m; disp('function Z = Gn(X)');...
  56.             disp('x = X(1); y = X(2);');...
  57.             disp('G = [2*x-4-y, 2*y-1-x];');...
  58.             disp('if norm(G)~=0,');...
  59.             disp('  Z = -G/norm(G);');...
  60.             disp('else');...
  61.             disp('  Z = -G;');...
  62.             disp('end;');...
  63. diary off;
  64.  
  65. % Remark. f.m Fn.m Gn.m grads.m are used for Algorithm 8.4
  66. % Remark. f(x,y) is used only to graph 2-dimensional functions.
  67. f(0,0); Fn([0 0]); Gn([0 0]);
  68. pause    % Press any key to continue.
  69.  
  70. clc;
  71. % The example is 2-dimensional, hence a surface z = f(x,y)
  72.  
  73. % can be graphed.  (Omit this section for higher dimensions.)
  74.  
  75. % The endpoints of the rectangle [a,b] x [c,d] are a,b,c,d.
  76.  
  77. % Place the maximum number of iterations in  max1   
  78.  
  79. % Place the tolerances in  delta  and  epsilon.
  80.  
  81. a  = 0;
  82. b  = 4;
  83. c  = 0;
  84. d  = 4;
  85. max1 = 50;
  86. delta = 1e-5;
  87. epsilon = 1e-7;
  88.  
  89. pause % Press any key to enter the starting vertices.
  90.  
  91. clc; clg;
  92. hs = (b-a)/16;
  93. ks = (d-c)/16;
  94. [Xs,Ys] = meshdom(a:hs:b, c:ks:d);
  95. Zs = f(Xs,Ys);  % User defined function.
  96. mesh(Zs);...
  97. title('To search for a minimum of z = f(x,y).');...
  98. shg; pause    % Press any key to continue.
  99.  
  100. clc;
  101.  
  102. % The gradient method is used to find the minimum
  103.  
  104. % of the n-dim function Fn(X). for convenience,
  105.  
  106. % For convenience, Fn(X) can be expressed using
  107.  
  108. % X = (x1,x2,x3,x4,x5,x6)  and the variables
  109.   
  110. % x = x1 , y = x2 , z = x3 , u = x4 , v = x5 , w = x6 
  111.  
  112. % You must supply the starting point P0(1:n)
  113.  
  114. P0 = [0.0   0.0];
  115.  
  116. pause    % Press any key to find the minimum of Fn(X).
  117.  
  118. clc; clg;
  119. a0  = 0;
  120. b0  = 3.5;
  121. c0  = 0;
  122. d0  = 2.5;
  123. hold off;
  124. axis([a0 b0 c0 d0]);...
  125. plot(P0(1),P0(2),'o');...
  126. hold on;...
  127. plot([a0 b0],[0 0],'b',[0 0],[c0 d0],'b');...
  128. xlabel('x');...
  129. ylabel('y');...
  130. title('Iteration points for the gradient method.');...
  131. grid;...
  132. shg; pause    % Press any key to continue.
  133.  
  134. [P0,y0,h,err,P] = grads('Fn','Gn',P0,max1,delta,epsilon,1);
  135.  
  136. format long;
  137. Mx1 = 'The coordinates of the point  P  are:';
  138. Mx2 = 'Error bound for the coordinates of  P:';
  139. Mx3 = 'The minimum value  F(P)  is:';
  140. Mx4 = 'Error bound for  F(P)  is:';
  141. clc,echo off,diary output,...
  142. disp(''),disp(Mx1),disp(P0),disp(Mx2),...
  143. disp(['± ',num2str(h),'      ± ',num2str(h)]),disp(''),...
  144. disp(Mx3),disp(y0),disp(Mx4),disp(['± ',num2str(err)]),...
  145. diary off,echo on
  146.